# 导入docplex 
from docplex.mp.model import Model
import numpy as np
# 创建模型对象 
model = Model() 
  #                8               12              16           20
PV1=[0,0,0,0,0,0,0,50,120,190,210,210,205,200,200,190,170,140,20,0,0,0,0,0] #1~24时刻光伏1发电
PV2=[0,0,0,0,0,0,0,0 ,0  ,100,150,200,220,220,200,200,180,150,50,0,0,0,0,0] #1~24时刻光伏2发电
PV3=[0,0,0,0,0,0,0,0 ,50 ,110,120,130,140,170,160,150,120,10 ,0 ,0,0,0,0,0] #1~24时刻光伏3发电
PV4=[0,0,0,0,0,0,0,0 ,0  ,20 ,50 ,60 ,90 ,100,100,100,110,115,50,0,0,0,0,0] #1~24时刻光伏4发电
PV5=[0,0,0,0,0,0,0,50,120,210,270,300,320,330,310,280,200,110,30,0,0,0,0,0] #1~24时刻光伏5发电
# c1=[1.05,1.03,1.00,1.05,1.05,1.10,1.05,1.00,0.85,0.86,0.50,0.48,0.43,0.55,0.45,0.53,0.60,0.80,1.00,1.00,1.00,1.02,1.02,1.01] #0~24时刻电网售电价格
c1=np.full(24,1.0)
c2=[1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.75,0.86,0.49,0.45,0.43,0.45,0.42,0.50,0.59,0.79,0.99,0.99,0.99,0.99,0.99,0.99] #0~24时刻ESP电价格
# c3=[0.45,0.45,0.45,0.45,0.45,0.46,0.47,0.49,0.50,0.46,0.48,0.45,0.45,0.45,0.45,0.45,0.45,0.45,0.48,0.49,0.50,0.45,0.45,0.45] #0~24时刻热网售热价格
c3=np.full(24,0.5)
c4=[0.46,0.46,0.46,0.46,0.46,0.47,0.49,0.50,0.50,0.47,0.49,0.46,0.46,0.46,0.46,0.46,0.46,0.47,0.50,0.50,0.50,0.46,0.46,0.47] #0~24时刻ESP热价格
# c5=[3.40,3.40,3.40,3.40,3.40,3.40,3.40,3.50,3.40,3.40,3.40,3.45,3.45,3.50,3.45,3.45,3.50,3.55,3.55,3.55,3.55,3.50,3.40,3.40] #0~24时刻气网售气价格
c5=np.full(24,3.5)
c6=[3.45,3.50,3.50,3.50,3.45,3.50,3.50,3.50,3.45,3.45,3.45,3.50,3.48,3.50,3.48,3.48,3.50,3.50,3.50,3.50,3.50,3.50,3.45,3.45] #0~24时刻ESP气价格
c_pv=0.039 #PV运行价格
a1=-0.5 #电力负荷相对于ESP售电价格的弹性系数
a2=-0.5 #电力负荷相对于ESP售热价格的弹性系数
a3=-0.5 #电力负荷相对于ESP售气价格的弹性系数
v1=2.5 #用户对电能偏爱参数
v2=2.1 #用户对热能偏爱参数
v3=2.3 #用户对气能偏爱参数
w1=0.0009 #用户对电能偏爱参数
w2=0.0011 #用户对热能偏爱参数
w3=0.0010 #用户对气能偏爱参数
#添加变量
P2=np.random.rand(0,10,24) #0~24时刻固定电负荷 
P3=10 #可调节电负荷上限
P4=[model.continuous_var(name=f"P4:{i}", lb=0, ub=P3) #43
             for i in range(24)] #0~24时刻可调电负荷
P5=np.array(PV1)+np.array(PV2)+np.array(PV3)+np.array(PV4)+np.array(PV5)
P1 = P2+P4-P5 #0~24时刻产消者电网负荷
P7=np.random.rand(0,10,24) #0~24时刻固定热负荷 
P8=5
P9 = [model.continuous_var(name=f"P9:{i}", lb=0, ub=P8) #46
             for i in range(24)] #0~24时刻可调热负荷
P6=P7+P9 #0~24时刻产消者热网负荷
P11=np.random.rand(0,10,24) #0~24时刻固定气负荷 
P12=5
P13 = [model.continuous_var(name=f"P13:{i}", lb=0, ub=P12) #49
             for i in range(24)] #0~24时刻可调气负荷
P10 = P11+P13 #0~24时刻产消者气网负荷
U1=sum(v1*x-(w1/2)*(x*x) for x in P1) 
U2=sum(v2*x-(w2/2)*(x*x) for x in P6)
U3=sum(v3*x-(w3/2)*(x*x) for x in P10)
U = U1+U2+U3
PV=sum(x*c_pv for x in P5)
I1=sum(c2*x for x in P1)
I2=sum(c4*x for x in P6)
I3=sum(c6*x for x in P10)

F=U-I1-I2-I3-PV
# 添加约束 
# for a,b,c,d in zip(P1, P2,P4,P5):
#   model.add_constraint(a-b-c+d==0) #42
for a,b,c in zip(P4,c1,c2):
    model.add_constraint(a-P3*(1+a1*(b-c)/b)==0) #44
# for a, b,c in zip(P6, P7,P9):
#   model.add_constraint(a-b-c==0) #45
for a,b,c in zip(P9,c3,c4):
  model.add_constraint(a-P8*(1+a2*(b-c)/b)==0) #47
# for a, b,c in zip(P10,P11,P13):
#   model.add_constraint(a-b-c==0) #48
for a,b,c in zip(P13,c5,c6):
  model.add_constraint(a-P12*(1+a3*(b-c)/b)==0) #50
# 添加目标函数 
model.maximize(F) 
# 求解优化问题 
# def cplex_main():  

#  冲突分析检查
# from docplex.mp.conflict_refiner import ConflictRefiner
# refiner = ConflictRefiner()
# res = refiner.refine_conflict(model)
# res.display()

solution = model.solve() 
# 获取结果 
if solution: 
    print(f"最优值为：{model.objective_value:.2f}") 
    print(f"P4的取值为:{P4}\n") 
    print(f"P9的取值为:{P9}\n")
    print(f"P13的取值为:{P13}\n")
else: 
    print("求解失败")


import matplotlib.pyplot as plt


A=[]
B=[]
C=[]
array = list(range(1, 25))
for a,b,c in zip(P4,P9,P13):
  A.append(a.solution_value)
  B.append(b.solution_value)
  C.append(c.solution_value)

plt.rcParams['font.sans-serif']=['Microsoft YaHei']  # 指定默认字体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号'-'显示为方块的问题

bar1=plt.bar(array,A, label='Electricity')
bar2=plt.bar(array,B,color='r',bottom=A, label='Thermal')
bar3=plt.bar(array,C,color='g',bottom=B, label='Gas')
plt.xticks(array)
# plt.grid(True)
plt.xlabel('Time/h')
plt.ylabel('Power/kW')
# plt.title('可调节能源策略')
plt.legend()
plt.show()
